Visualizing properties of population dynamics in Julia

by Rene L. Principe Jr.

PROJECT OBJECTIVES

1. To use Julia to plot the time development of a growing population as shown in Fig. 3.1 using the iteration Eq. 3.1 introduced by P.F. Verhulst [1].

$x_{n+1} = 4r(x_n)(1-x_n) = f(x_n), \quad n = 0, 1, 2, 3, ..., $

2. To replicate the following plots from [1]

r ≃ 0.934

r = 0.874640 (superstable four-cycle)

r = 0.9642 (chaotic three-cycle).

PROJECT LEARNINGS

This project mainly demonstrates the properties of the elementary iteration process of a nonlinear function that is the logistic map. With proper visualizations, we can show how sensitive chaotic systems are to their initial conditions [G3: Implement known numerical methods/algorithm for solving or simulating basic physical models using Julia as programming language]. Since this is an iterative process, essential steps shall be benchmarked and the runtimes will be compared to the implementation carried out in Python [G1: Identify basic principles in high-performance computation especially those accessible via Julia as a programming language].

References

[1] W. Kinzel and G. Reents, Physics by Computers: Programming Physical Problems using Matematica and C (Springer-Verlag, Berlin Heidelberg, 1998). URL: https://doi.org/10.1007/978-3-642-46839-1

[2] https://en.wikipedia.org/wiki/Logistic_map

Deterministic chaos is defined as the unpredictable behavior of a well-defined function given a sensitive initial information. This is carried out through the process of iteration, or essentially repeating the same operation unto itself - in this case, the quadratic function. On the visualizations in this project, we will see how regular behaviour turns chaotic given a quantitative value [1].

In this population dynamics report, an equation introduced by P.F. Verhulst shall be used which is given by the iteration

$x_{n+1} = 4r(x_n)(1-x_n) = f(x_n), \quad n = 0, 1, 2, 3, ..., $

Time Development of a growing population

To model a simple time development of a growing population, the growth of a population density is described by $4rx_n$. In real life, food supply is actually limited, hence we must incorporate the nonlinear effect that is $-4rx^2_n$. Defining the function $f(r,x)$, we have

where the solution is an approximation of a differental equation

$x(n) = \frac{4r-1}{4r + \text{const·exp}((1-4r)n)}.$

Now, we setup the code to show us the time development and display its plot. The function takes in $N$ or the number of iterations desired and the parameter $r$, which will be varied to explore different scenarios.

r < 1/4 → x∞ = 0

To visualize the behavior of the equation, an animation was made for varying r values. Here we can see various transition points.

r > 1/4 → x∞ = 1 - 1/(4r)

Evidently, the phase transition is at $r_o = 1/4$ where below this value, the attractor approaching zero implies vanishing population and above it, the population converges to a constant, that is $x∞ = 1 - 1/(4r)$.

Interestingly, for larger values of $r$, we can see that there are multiple attractors, meaning that the population coverges to multiple values which are known as stable fixed points. For $r = 0.87$, the $x_n$ values moves toward an attractor of period 4. Also, we can see how the $x_n$ values jumps back and forth but the pattern is starting to become very evident. It might seem like random but in fact, these points follow a well defined parabola. This behavior is what comprises deterministic chaos.

The perturbation essentially explodes at $r > 3/4$ and the two-cycle eventualy becomes four-cycle. As it reaches $r_o = 0.89$, the system becomes way to chaotic, which implies that the doubling of periods has reached infinity.

The function and its two- and fourfold iterations

It was shown previously how at $r = 0.87$, $x_n$ values tended towards the attractor of period 4. This implies that the function has been iterated for four times already (fourfold iteration, $f4(x)$). The plot above shows how f4(x) intersects the $f0(x)$ but only four are stable points.

The Logistic Map

The following plots are logistic maps containing 1000 $x_n$ values from the iteration equation, exlcuding the first 100 points. These corresponds to 1000 vertical lines for each r value, which was also varied from 0 - 1.

The plot below highlights how attractors are properly visible at some $r$ values, and then some "forking" of these stable orbits begins as attractors converge to 2, and then to 4, 8, etc. until the forking becomes chaotic. It best highlights the concept of chaos because treating parameter $r$ as our initial condition, we have shown how it yielded wildly differing results, hence, it is very sensitive to the initial condition.

Looking closer, we can clearly see where the stable period transitions to the chaotic region.

Upon a much closer inspection, stable regions are in fact present within the chaotic region.

Histograms

Another way to visualize the property of chaos is by looking at the frequency distribution of the $x_n$ values where local peaks would imply the approximate period.

Hence, for $r = 0.934$, the approximate period if 5, referring to the five peaks in the histogram.

Cobweb plot | Verhulst diagram

In here, we construct the cobweb plots. The animation displayed below shows how the algorithm works. It is able to show if a point is stable whenever it spirals inward. Otherwise, unstable fixed points would yield a cobweb plot where the spirals are oriented outwards.

For example, shown above is the cobweb plot at fixed point $r = 0.934$. As per the bifurcation diagram and the histogram, we know that it behaves chaotically. Its manifestation in this cobweb plot is through the complex closed loops exhibited by the animation.

To characterize the loop behavior of the stable and unstable fixed points, we animate the cobweb plot for $r$ values from 0.5 - 1.

Shown above is the transition from how the closed loop spirals inwards from $r = 0.5-0.8$, it became super stable at $r ~ 0.8$, and it quickly transitioned into a chaotic period as the loops spiraled outwards and it becomes more complex.

The code below returns a cobweb plot after 100 iterations, given an $r$ value and a starting point. We can see how it

Combining the bifurcation plot, histogram, and the cobweb plots, we examine three fixed points.

r = 0.874640 [Superstable four-cycle]

r = 0.934 [Ultra chaotic five-cycle orbit]

r = 0.9642 [Three chaotic bands near the three-cycle stable fixed point]

In summary, the properties of population dynamics were successfully visualized through the time evolution plots, logistic map, histogram, and the cobweb plots. Animations were also implemented to highlight the dynamic behaviour the this population model. Hence, the third goal [G3: Implement known numerical methods/algorithm for solving or simulating basic physical models using Julia as programming language] for this course was succesfully adopted in this project. In comparison to Python, producing animated plots was a lot easier in Julia. Due to time and resources (PC) constraint, I wasn't able compare the benchmark scores in Python and Julia but, @benchmark was carried out on some parts of the code which shows how the algorithm performs and works.